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1. Introduction 

Spectral methods consist of expanding the solution to a problem in terms of basis functions which are 
global, infinitely differentiable and preferably orthogonal [1,2]. This choice of basis functions is what distin- 
guishes them from the finite difference and finite element methods. In the case of finite element methods, the 
domain is divided into small elements, and a basis function is specified in each element They are thus local 
in character. The case with finite difference methods is similar. 

In addition to the basis functions, a key element of spectral methods is the set of test functions or 
weight functions. The test functions are used to enforce minimization of the residual resulting from the sub- 
stitution of the series expansion of the solution into the differential equation. The choice of test functions 
distinguishes between essentially two types of spectral methods — spectral Galerkin, and spectral collocation. 
In the Galerkin approach, the test functions are usually the same as the basis functions, w'hereas in the collo- 
cation approach the test functions are translated Dirac delta functions. In other words, the Galerkin approach 
satisfies the differential equations in the least square sense. In the spectral collocation approach the equations 
are satisfied exactly at the selected, so-called collocation points. It should be noted that the basis functions are 
employed solely for the purpose of approximating derivatives. This approach, which became feasible with the 
advent of computers, is the easiest and the most efficient for nonlinear problems, . and is the focus of the 
present discussion. 

2. Basic Aspects 

The principal advantage of spectral methods lies in their potential for rapidly convergent approxima- 
tions. In practical terms, it means that they achieve accurate results with substantially few-er points than are 
required by typical finite difference methods. Suppose u A , is a numerical approximation to a function u(x). 
With a given set of basis functions <£„, it takes the form 

N 

udx) = X «« (1) 

The expansion coefficients a r are obtained by enforcing the condition 

u N (xj) = u(xj). (2) 

where xj are the selected, so-called collocation points, which are usually the extrema of Figure 1 pro- 
vides a graphic distinction between a second-order accurate central difference approximation and a Legendre 
spectral approximation to the first derivative of the function 


on [-1, 1] 


u(x) = 1 + sin(2:u + —) 
4 


whose values are given at a finite number of grid points. The finite difference approximation for the deriva- 
tive at the origin, for instance, is estimated by interpolating a parabola through the origin and the two adja- 
cent points, and is thus local in character. The spectral approximation estimates the derivative of the original 
function by the derivative of the polynomial which interpolates all the available points. Note that the error of 
the finite difference discretization decreases as 1 IN 2 , whereas the error of the spectral discretization decreases 
exponentially. In the case of a differential equation, a further step is involved, that of finding an approxima- 
tion for the differential operator in terms of the grid point values 

Another advantage of spectral methods is their minimal phase error. Consider the periodic solution to 
the problem u t + u x = 0 with w(x,0) = sin(;tcos(x)). Figure 2 shows the lagging phase of the finite difference 
solution, while the Fourier spectral solution has zero phase error. A fourth-order Runge-Kutra method is used 
for temporal discretization in both the cases. For realistic problems with variable coefficients or nonlinear 
terms, the phase error for spectral methods is, of course, nonzero, but still relatively small. 

These are some of the essential aspects of spectral methods which make them the prevailing tool in the 
study of stability, transition, and turbulence. Some of the drawbacks which have inhibited their wider use are: 
1) time-step restriction imposed by the standard spectral grid, 2) sensitivity to singularities and 3) restriction 
to simple geometry. Progress has been made on all counts [3]. The present work will be confined to the 
recent developments in overcoming the first obstacle. 


3. Iterative Spectral Methods 

For evolution problems, explicit time-stepping can be extremely inefficient. This is so because the typi- 
cal time-step limitation for spectral methods is proportional to 1 IN 2 for the advection equation and l/N* for 
the diffusion equation (where N is the number of modes) [1]. Hence, implicit time-stepping becomes a neces- 
sity. This results in a set of algebraic equations, which are in general amenable to iterative solution tech- 
niques only. Also, elliptic equations governing practical problems virtually require implicit iterative tech- 
niques. As the condition number of the relevant matrices are large, preconditioned iterative schemes including 
multigrid procedures are the attractive choices. In this section, we discuss the fundamentals of iterative spec- 
tral methods with reference to an elementary example, and then their application to the three dimensional 
incompressible Navier Stokes equations for the study of stability and transition in a channel and a boundary 
layer. 

For the purpose of illustration let us consider the equation. 


periodic on [0,2 ti]. For the Fourier method, the standard choice of collocation points is 

_ 2tc/ 


* m 'N- 


j = 0, 1,2 N- 1, 


( 3 ) 

( 4 ) 


Setting uj = u{xj), the discrete Fourier series for u may be represented by the discrete transform pair 
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The expression for the derivative u x at the collocation points is 


« A ) = Z <> v^' . 

tf . 



j - 0, 1, .... N-l 


Thus the Fourier collocation discretization of the equation may be written 


LU = (5) 

where [/ = (hq, uj, w^), F = C/o»/ii ^-i)» and £ = C~ l DC with C being the discrete Fourier transform 
operator, CT 1 the inverse transform, and D the diagonal matrix denoting the first derivative operator in the 
Fourier space. Specifically, 

(HW2) 

_ jfj j, k - 0, 1 N-l (6) 

^jk "" c » 

and 


Djj = i (j — N/2) 


for j - 1, 2 N-l (7) 


= 0 for j - 0 

Hie eigenvalues of L are Up) = ip, p = - N/2 + 1, .... N/2 — 1, and the largest one grows as N/2. A 
preconditioned Richardson iterative procedure for solving Eq. (5) is 

V <- V + to/T 1 (F - LV) (8) 

where the preconditioning matrix H is an approximation to L, is sparse, and is readily invertible. An obvious 
choice for H is a finite difference approximation Lpp to the first derivative. With the various possibilities for 
L F q, the eigenvalue spectrum of L FD L is given in Table 1. Apparently, the staggered grid leads to the most 
effective treatment of the first derivative. This kind of preconditioning was successfully used in the semi- 
implicit time-stepping algorithm for the Navier Stokes equations discussed in the section on Navier Stokes 
Algorithms. The eigenvalue trends of that complicated set of vector equations are surprisingly well predicted 
by this extremely simple scalar periodic problem. 

Next, let us consider the second order equation 

on [0,2 ti] (9) 

with periodic boundary conditions. A Fourier collocation discretization of this equation is the same as Eq. (5) 
except for the diagonal matrix D which represents now the second derivative operator in the Fourier space. 
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;= i,2. 


A— 1 


= 0, 


j = 0 


The eigenvalues of L are \{p) = /t 2 , /> = -Nil +1, .... A72-1. To make the case for the multigrid procedure 
(consisting of a fine-grid operator and a co arse-grid correction) as a preconditioner, we assume H to be the 
identity matrix I in the iterative scheme (8). The iterative scheme is convergent if the eigenvalues, (1 - (oX), 
of the iteration matrix [/-coL] satisfy 

|1 - coX| < 1. 

Each iteration damps the error component corresponding to X by a factor v(X) = |1-G)X|. The optima] choice 
of X is that which balances damping of the lowest-frequency and the highest-frequency errors, i.e., 

(1 - ttXmix) = - 0 - ttXmin) 


This yields 


= 


(^max ^min) 


and the spectral radius 


Rsg = 


C^-max ^ttud) 
(^4nax ^ min ) 


In the present instance, X^ = N 2 /4 t X^ = 1, and thus = 1 - 8 IN 2 . This implies order N 2 iterations are 
needed for convergence. This poor performance is due to balancing the damping of the lowest frequency 
eigenfunction with the highest-frequency one. The multigrid procedure exploits the fact that the lowest- 
frequency modes (|p| < A74) can be damped efficiently on coarser grids, and settles for a relaxation parameter 
value which balances the damping of the mid-frequency mode (^J = Nl 4) with the highest-frequency one 
(lp| = Nil). Table 2 provides the comparison of single-grid and multigrid damping factors for AM>4. The 
high frequencies from 16 to 32 arc damped effectively in the multi grid procedure, whereas the frequencies 
lower than 16 arc hardly damped at all. But then some of these low frequencies (from 8 to 16) can be 
efficiently damped on the coarser grid with A r =32. Further coarser grids can be employed till relaxation 
becomes so cheap that all the remaining modes can be damped. In concrete terms, the ingredients of a mul- 
tigrid technique are a fine-grid operator, a relaxation scheme, a restriction operator which interpolates a func- 
tion from the fine grid to the coarse grid, a co arse-grid operator, and a prolongation operator interpolating a 
function from the coarse grid to the fine grid. The fine grid problem for the present example may be written 

L f U f = F f (10) 


Let denote the fine-grid approximation. After the high-frequency content of the error V^— (/has been 
sufficiently damped, attention shifts to the coarse grid. The co arse-grid problem is 


1 + 



L e ir = F* 


(ii) 


where 

/* = * [Ff-lJvf], 

R being the restriction operator. After a satisfactory approximation V* is obtained, the coarse-grid correction 
(V* — RV*) is interpolated onto the fine grid by the prolongation operator P, yielding the corrected fine-grid 
solution 

vf*~ V'+PiV'-RV') (12) 

The details of spectral multigrid techniques are furnished in [4]. Spectral multigrid techniques have been 
used to solve a variety of problems including the transonic full potential equation [5,6]. Additional applica- 
tions of spectral methods to compressible flows are described in [7]. In the next section, we describe a mul- 
tigrid algorithm for the incompressible Navier Stokes equations. 

4. Navier Stokes Algorithms 

This section is devoted to a description of recently developed algorithms for the simulation of instability’ 
and transition to turbulence in a flat-plate boundary layer. These algorithms deal with the primitive variable 
formulation of the Navier Stokes equations, and are based on the iterative methods discussed above in the 
simplest context They are capable of handling geometric terms and variable viscosity. 

The Navier Stokes equations in the so-called rotation form are 

q t = q x co + V * QiV?) — VP in Q 

V • q = 0 in Q 


and 


q{Xf 0) = <7oto 


in Cl (13) 


q = g on dCl 

where q = (w,v,w) is the velocity vector, to = V x q the voracity, P = p + 1/2 |< 7 i 2 the total pressure, p the 
variable viscosity, Cl the interior of the domain, and dCl its boundary. In the stability and transition problems 
under study, the domain Cl is cartesian and semi-infinite: periodic in the two horizontal directions (x,z), and 
bounded by a wall at y= 0. Fourier collocation can be used in the periodic directions (x,z), and Chebyshev 
collocation is used in the vertical (y) direction. The collocation points in the periodic directions are given by 
a relation similar to Eq. (4). The vertical extent of the domain 0 < y < «© is mapped onto — 1 < £ < +1. The 
velocities are defined and the momentum equations enforced at the points 

Sy= cos( -j^), j-0,1 Ny 

The pressure is defined at the half points 



= cos[ 


n (ffl/2) 
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j - 0, l N y * l 


5 



and the continuity equation is enforced at these points. The staggered grid avoids artificial pressure boundary 
conditions, and precludes spurious pressure modes. 

After a Fourier transform in x and z, the temporal discretization (backward Euler for pressure, Crank- 
Nicolson for normal diffusion, and third or fourth-order Runge-Kutta for the remaining terms) of Eqs. (13) 
leads to 

[/ - MDM) Q + At A 0 Vn = Q t (14) 


-A+VQ = 0 

where 

Q = {9o +I . $*' 

n = {PV2> Pin PT\ai (15) 


V 



M is the Chebyshev derivative operator, D the diagonal matrix with l/2pAr as its elements, and Aq is the 
interpolation operator from the half points to cell faces, A+ vice versa. Obviously, the equations for each pair 
of horizontal wave number (k^kj are independent, and they can be written as the system 


LX = F 


where X = [Q, II] . The iterative solution of this equation is carried out by preconditioning the system with a 
finite difference approximation on the CKebyshev grid, and applying a standard iterative technique such as 
Richardson, minimum residual or multigrid [8J. 

The method described above solves the implicit equations together as a set The extension of this 
method to the more general cases of interest such as those involving two or more inhomogeneous directions 
is not straightforward. An alternative is the operator-splitting technique or the fractional step scheme [9]. This 
method yields implicit matrices which are positive definite and are easily amenable to iterative methods. In 
the first step, one solves the advection-diffusion equation 

q* = < 7 * X CD* + V * (pV^*) (1 6) 


subject to the initial and boundary conditions 


q* ~ g* on dft 

Note that g* has yet to be defined. In die second step, one solves for the pressure correction 

q = VP** 


( 17 ) 


V • q** = 0 
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subject to the conditions 


q**(xS) = q in n 
q ** • n = g • n in 

where n is the unit normal to the boundary. Further, the tangential component of the Eq. (17) holds on the 
boundary, i.e., 

q** • t = - VP ** . t in an 

where t is a unit tangent vector to the boundary. Now g* is defined [9] as (using Taylor expansion in t) 

g* * n = (g* + Ar g*) * h 

g** + 

Eq. (16) is discretized in the usual spectral collocauon manner. After a temporal and spatial discretization of 
Eq. (17), the boundary conditions are built into the relevant matrix operators, and then a discrete divergence 
is taken. This results in a discrete Poisson equation (with as many algebraic equations as unknowns) for pres- 
sure, which can be solved by standard iterative techniques including the multigrid method. 

5. Applications 

These algorithms have been used to study the incipient stages of the transition process in channel flows 
[8] and parallel boundary layer flows [10]. Some representadve results are provided here. The channel flow 
results pertain to the secondary instability associated with the so-called center modes. Unlike the Tollmien- 
Schlichting modes (sometimes alluded to as wall modes), the center modes always decay with a rather high 
decay rate. Their phase velocity is near unity and their maxima occurs near the center of the channel. The 
simulation had a Reynolds number of 5000 based on half-channel width, and the initial conditions consisted 
of a 20% amplitude two-dimensional center mode with two 1.5% amplitude skewed modes. The harmonic 
contents of the solution were monitored, and die grid was refined as deemed necessary. The finest grid was 
144x96x108. Plotted in Figure 3 are stream wise vcrticity (left side) and span wise vorticity (right side) con- 
tours on the streamwise planes at x * 3/8, 7/16, 1/2 and 9/16 of the fundamental wavelength. The so-called 
peak plane would intersect these streamwise planes along a vertical line in the center of the frames. The 
structure of the vortex loop can be deduced from these plots. It differs in detail from that of the wall modes. 
The vortex structures are significant over only a small portion of the wavelength in the streamwise direction, 
whereas in the case of wall modes they cover almost the whole wavelength. Furthermore, the pinching of the 
vortex loop in the peak plane appears to be less acute in the case of the center modes. The harmonic history 
is displayed in Figure 4. The evolution of the secondary instability is apparent and it appears similar to that 
of the Tollmien-Schlichting modes. What is more interesting is the steep growth of the (0,2) and (2,2) modes 
which may lead to a strong tertiary instability. To resolve it in detail would require an even finer grid. 

The parallel heated water boundary layer cases had a Reynold number of 1100 based on the displace- 
ment thickness, and the initial amplitudes of the two-dimensional and three-dimensional Tollmien-Schlichting 
waves were 2.7% and 0.4% respectively. 
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Three different situations were srudied: 1) uncontrolled case, 2) heated fixed temperature case, and 3) 
heated active temperature case. In the heated fixed temperature case, the temperature was kept fixed at the 
initial value pertinent to the mean flow conditions, and the temperature evolution was totally neglected. In the 
heated active temperature case, the temperature evolution was taken into account by solving the temperature 
equation along with the momentum equations. In both the cases the wall temperature was 2.75% above the 
free stream temperature. Figures 5,6 and 7 display the harmonic histories. The fixed temperature case over- 
predicts the weakening effect of heating on the secondary instability. Figure 8 shows the spanwise voracity 
contours on the peak plane. In the uncontrolled case (Figure 8 top left) a kink develops in the high-shear 
layer at time / equal to three Tollmien-Schlichting periods. It is generally accepted that a irrevocably quick 
succession of events follows thereafter leading to a turbulent spot formation. Heating the wall to 2.75% 
above the free stream temperature diffuses the high- shear layer as is obvious from Figure 8 (bottom left). 
However, within the subsequent one and one fourth period, turbulent spot formation appears to become 
imminent (Figure 8 top right). In the fixed temperature case, it is clear from Figure 8 (bottom right) that the 
high-shear layer formation is mellowed down even up to four and one fourth periods. This shows that the 
effect of temperature evolution is significant and deleterious in the nonlinear regime whereas it is quite negli- 
gible in the linear regime. 
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Table 1. Preconditioned Eigenvalues for One-dimensional First 
Derivative Model Problem 


Preconditioning 

Eigenvalues 

Central Differences 

kAx 

sin(JkAx) 


■ ,■ \m * (2jc/3) 

sm(JtAx) 

High Mode Cutoff 

0 (2n/3) < |*Ax| <S it 

One-sided Differences 

-tikAxfl) kAx/2 

sin((JkAx)/2) 

Staggered Grid 

(kAx)!2 

sin((JfcAx)/2 


Table 2. Damping Factors for N - 64 


p 

Single-Grid 

Multigrid 

i 

.9980 

.9984 

2 

.9922 

.9938 

4 

.9688 

.9750 

8 

.8751 

.9000 

12 

.7190 

.7750 

16 

.5005 

.6000 

20 

.2195 

.3750 

24 

.1239 

.1000 

28 

.5298 

.2250 

32 

.9980 

.6000 
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Figure 2: Finite difference (FD2) and Fourier spectral (FS) approximations 
after one period to a simple wave equation whose exact solution is represented 
by the curve. 
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Figure 5: Harmonic history for a Re = 1100 boundary layer. 
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Figure 7: Same as Figure 6 but for fixed temperature. 
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